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A new theoretical approach is described for the inverse self-assembly problem, i.e., the reconstruction of the 
interparticle interaction from a given structure. This theory is based on the variational principle for the 
functional that is constructed from a free energy functional in combination with Percus’s approach [J. Percus, 
Phys. Rev. Lett. 8, 462 (1962)]. In this theory, the interparticle interaction potential for the given structure 
is obtained as the function that maximizes the functional. As test cases, the interparticle potentials for 
two-dimensional crystals, such as square, honeycomb, and kagome lattices, are predicted by this theory. The 
formation of each target lattice from an initial random particle conhguration in Monte Carlo simulations with 
the predicted interparticle interaction indicates that the theory is successfully applied to the test cases. 


I. INTRODUCTION 

In condensed phases, diverse interactions between the 
constituents (such as atoms, molecules, micelles, nano- 
and microparticles) produce a variety of simple and 
complex structures. Besides a plethora of the crystal, 
liquid crystal, and quasicrystal structures formed by con¬ 
ventional atoms and molecules, square and honeycomb 
lattices of colloidal nano crystals,^ kagome lattice of tri¬ 
block Janus particles,^ and quasicrystals of dendrimers^ 
and binary nanocrystals^ are additional experimental 
examples of non-trivial structures. Molecular simulations 
for systems with several interparticle interactions have 
been performed to enumerate the stable structures; it has 
been found that particles with a rather simple interaction 
potential can assemble into complex structures, e.g., 
hard sphere (HS) plus linear ramp,® HS plus square 
well,® HS plus square shoulder,^ and Lennard-Jones 
plus Gaussian® (LJG) potential particles assemble into 
a two-dimensional (2D) quasicrystal. Density functional 
theory (DFT), another tool to hnd the thermodynami¬ 
cally stable phases, has shown that in three dimensions 
(3D), the LJG fluid forms a body centered cubic (bcc) 
crystal as a stable phase.® 

When a structure of a condensed phase is known from 
experiments or is designed artificially, determination of 
the interparticle interaction is more efficient by using the 
inverse approach, in which the interparticle interaction 
is derived from the given structure, than using the 
forward approach of exhaustive search by the molecular 
simulations or DFT. The inverse statistical mechanical 
method,^® which is the most elaborated inverse approach, 
has been successfully used to design the interparticle 
interaction potentials that generate the square, hon¬ 
eycomb,kagome, and rectangular^®d4 lattices in 
2D, and the simple cubic, bcc, simple hexagonal, 
diamond,17 wurtzite lattices^"^ in 3D. In the 

inverse statistical mechanical method, the interparticle 
interaction is optimized with respect to the energy and 
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mechanical stability criteria via molecular simulations. 
The studies on the inverse problem at finite temperature 
have mostly been based on the molecular simulations 
at finite temperature. (It was suggested in Ref. 18 
that the molecular simulations at finite temperature may 
not be necessary. This is based on the observation 
that the target structure shows good stability at finite 
temperature if the interaction potential is optimized 
at zero temperature via the minimization of a specific 
simulated annealing energy defined in Ref. 16. However, 
this strategy is also dependent on computer simulations.) 
For the liquid target phases, the reverse Monte Garlo 
(MG) method,^®’®® which also uses the MG simulation in 
the optimization, is another successful inverse approach. 


Here the aim is to formulate a theory to reconstruct the 
interparticle potential that can stabilize a given target 
structure (more specifically, the single- and two-particle 
distributions) at finite temperature. In this paper, I 
present a new simulation-free inverse method that is a 
variational method based on the free-energy functional 
theory®^ akin to DFT. The interparticle interaction func¬ 
tion is defined as the function that gives the maximum 
of the functional. I applied this method to the square, 
honeycomb, and kagome lattices, and obtained the corre¬ 
sponding interparticle potentials. The potentials are then 
used in a series of simulated annealing MG calculations 
starting from random configurations. In most cases, the 
resulting solid contains a few grain boundaries and many 
defects. However, it is found that for each predicted 
potential, in at least a few percent of the simulations 
the particles spontaneously form the appropriate target 
lattice with a small number of defects. Although the 
small success rate in the simulations indicates that the 
interaction potentials obtained here are not entirely opti¬ 
mized as the potentials obtained by the inverse statistical 
mechanical method, the observed self-assembly into the 
target lattice implies that the method introduced here 
is another promising approach to the inverse problem of 
the self-assembly. 
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II. THEORY 


We consider a single component system in a d- 
dimensional space. The system comprises particles in¬ 
teracting with a pairwise-additive potential v(x). The 
grand potential of the system, with the temperature T, 
chemical potential /i, and volume V, iir the presence of 
an external field is defined as 

( 1 ) 


where /? = l/feeT is the inverse temperature, ^p{x) = 
— (jjextix) is the intrinsic chemical potential, and S[(/5] 
is the grand partition function 
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where and A denote the coordinate of the ith particle 
and the de Broglie thermal wavelength, respectively. 

The key functional for this approach is 


A[p,^] = + 


p{x, [>f])tlj{x)dx, 


( 3 ) 


where p{x^ [</?]) denotes the single-particle density in the 
presence of the intrinsic chemical potential p>{x) and 
ipi^x) is an independent function. The maximum of 
A[yO,'(/'] with respect to ip, i.e., 

A[p\= sup A[p,iP], (4) 

p 


is the intrinsic free energy, the Legendre transformation 
of the grand potential Al['?/’]• The maximum of A[p,il)] is 
achieved when ip is equivalent to the intrinsic chemical 
potential p, i.e., A[p] = A[p,tp\. The inequality relation 

A[p, ip] > A[p, Ip] for any ip{x) (5) 

provides us the variational method^^ to determine the 
intrinsic chemical potential p{x) that gives rise to a 
given density profile p{x, [(/j]): the ip{x) that maximizes 
the functional A[p,ip] is the intrinsic chemical potential 
p{x). The variational method for A[p,ip], in which the 
proper p(x) is obtained for a given p{x), can be viewed 
as the inverse of the DFT,^^’^^ in which the proper p{x) 
is obtained for a given p{x). 

The variational method for A[p, ip] can be adapted to 
obtain the information about the interaction potential 
v{x) using Percus’s idea.^^’^^ If a particle is fixed at the 
origin of the coordinate system, the remaining particles 
feel the external field v{x) due to the fixed particle in 
addition to the original external field (pex.t{x). In this 
situation, the single-particle density becomes p{x, 
where pf,,^{x) = p — [(^ext(ic) -l-u(a:)] = p{x) — v{x) is the 


intrinsic chemical potential in the presence of the fixed 
particle. Percus^^ showed that p{x, is related to 

the single- and two-particle density in the absence of the 
fixed particle by 




p(o, M) 


( 6 ) 


where p^‘^\x,x', [(^]) is the two-particle density between 
X and x' in the absence of the fixed particle. As discussed 
above, the function ip that maximizes A[p[x, [pf,^),ip] is 
PfLx{x). Combining the variational method for A[p,ip] 
with Percus’s idea, we find that the function ip{x) that 
maximizes A\^p^‘Ajp^ip] is Thus, for a given set 

of p{x) and p^‘^\x, x'), we can obtain the interaction po¬ 
tential v{x) = p{x) — pfixix) through the maximization 
of A[p^'^'>/ p, Ip] with respect to ip{x). 

The inhomogeneous version of Percus’s relation (6) 
can be derived along the same line as the derivation of 
the homogeneous version in Ref. 22 by simply removing 
the assumption of the homogeneity. In an inhomoge¬ 
neous and symmetry-broken phase, however, the naive 
definition of the distribution functions loses its unique¬ 
ness. In some systems, for example, a ferromagnet, 
the uniqueness of the order parameter is recovered by 
applying a field that breaks the symmetry of the original 
system, and the order parameter without the field can be 
obtained by taking zero-limit of the symmetry-breaking 
field followed by the thermodynamic limit. I do not know 
whether the same procedure applies to the distribution 
functions of inhomogeneous fluids, thus cannot provide 
a rigorous proof of Percus’s relation for inhomogeneous 
cases; however, I postulate the validity of Eq. (6) in this 
paper. 


III. APPLICATIONS TO 2D CRYSTALS 

The free-energy functional method introduced in Sec. 
II is applicable to any 2D and 3D system that is 
characterized by its single- and two-particle densities. 
In this paper, as test cases, I use 2D crystals such 
as square, honeycomb, and kagome lattices as target 
structures. While the interaction potentials that stabilize 
these target crystals have already been found in previous 
studies^dO-14 gj.g currently of little novelty; never¬ 
theless, these crystals still serve as good test cases. 

From the variational principle for A^p^^'^/p,ip], it 
follows that the intrinsic chemical potential pfix^(x) 
is the solution of the Euler-Lagrange equation 
6A[p^‘^'>/P,ip]/Sip{x) = 0, i.e., 

p^'^\x,x') (50 [' 0 ] 

p(x') 5ip{x) 

A direct solution of this equation is computationally 
demanding. Therefore, in this paper, I use a trial 
function r’tr(35, {a^i}) with a set of variational parameters 
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{xi} as the interaction potential. The function ^/’(a;) 
then becomes -iptrix, {xi}) = fi — Vtr{x, {xi}), and the 
functional /p, ip] is reduced to the function of {xt}. 

The set of parameters {xi} that maximizes A j p, iptP] 
defines the required interaction potential. 

In the absence of the external field, the equilibrium 
configuration is determined solely by Pv(x), because 
the grand partition function (2) depends on v{x) only 
through f5v{x). Thus, the optimal interparticle potential 
is linearly dependent on the temperature. 

The grand potential 17 can be expanded in a 
functional Taylor expansion. Since iptrix) for the typ¬ 
ical interparticle interactions diverges around x = 0 
due to the short range repulsion, the activity z{x) = 
exp[/3'!/)tr(a:)]/A^ is easier to handle than iptr{x) itself. 
The functional Taylor expansion of the grand potential 
in powers of the deviation of activity Az(x) = z(x) — zq 
is 


pnllptr] = - ^ 


n—0 


(5” ln5[(p] 


Sz{xi)Sz(x2). ■. Sz{xn) 


z=zo 


'^Az{xi)dx^, 

( 8 ) 


where zq = j}? is the activity in the absence of 

the fixed particle. The functional derivatives of InS 
with respect to zix') can be expressed in terms of the 
multiparticle distribution functions.Substituting (6) 
and (8) into (3), we obtain 

PA[p^'^'> / P,ipt^^ =PQ,Yp\+P j ^ ^|^j^V tr(a^)da; 

- J p{x)Cix)dx-^ jj [p^‘^\x,x') 

- p{x)p(x')\C,{x)C,{x')dxdx', (9) 

up to the second order in Az{x), where C(®) = 

Az{x)/zo = exp[—/3utr(a;)] — 1- The approximation 
of the truncated Taylor series is justified if C,{x) is 
small. Therefore, when Vtrix) has a negative value, this 
approximation may break down at very low temperature. 
Hereafter, I restrict the consideration to a case in which 
the minimum of fdvtrix) is —1. 

As the trial function vt^ix), I use the LJG potential, 
which is an isotropic potential constructed by adding a 
Gaussian function to the Lennard-Jones potential.® Here, 
the trial function is the LJG potential with positive and 
negative Gaussians: 


■Ctr(£c) = 



of 1 V 


(ax — Xi)"^ 

^[axj 

— Cl exp 

[ 2a2 J 

(ax — X 2 )^ 
2al 


(10) 


where Ci and C 2 denote the Gaussian functions’ width, 
depth, and height, respectively. In this paper, I use the 


parameters (T„ = ^0.02, Ci = 1, and C 2 = 4. The 
variational parameters Xi {i = 1,2) are used to adjust 
the positions of the Gaussian functions. The values 
of the remaining parameters a and e are determined 
numerically so that the first minimum of Utr(a:) will be at 
the nearest neighbor atomic distance and the value of the 
global minimum will be —1. Gonsidering these criteria for 
the trial function, the functional A[p^'^'> jp, iptr] is reduced 
to a function A (a; 1 ,^ 2 ). 

The atomic positions in the target perfect crystals, 
{a^}, are the input for the free-energy functional method. 
I use units for which the nearest neighbor atomic distance 
is unity. In the present paper, I assume that the ith 
particle fluctuates around and is independent of the 
other particles, and that its fluctuation is given by 
a Gaussian distribution. The single- and two-particle 
distribution functions are then 

Pix) =Y^fiix-a,), (11) 

i 

p^^\x,x') =Y^fi{x-a,)^fi{x'-aj), (12) 

i 

where fi{x) denotes the Gaussian distribution 

fi{x) = ^exp 

The resulting interaction potential depends on the stan¬ 
dard deviation a. In principle, any value of a is 
permissible in this method. However, because p{x) and 
p^^'>{x,x') with sharp peaks are not preferable for the 
accuracy of the numerical integration in (9), the broadest 
but physically acceptable distribution fi{x) is better. 
Therefore, in this paper, cr is set to 0.15, which is the 
typical value of the Lindemann ratio at crystallization. 
Since p^^\x, 0)/p{0) remains finite near a; = 0, where 
ip^x{x) diverges rapidly, the approximations (11) and 
(12) cause the second term at the right-hand side of (9) 
to diverge. This is caused by neglecting the correlation 
between the particles in the derivation of (12). To 
include the non-negligible particle correlation due to 
the strong repulsion between the particles separated by 
short distances, I modified (12) by multiplying it by 
exp[-/3(l/|a;|i2_l)] for] a:| < 1. The modification factor 
is proportional to the ideal gas density in a repulsive 
external field l/|a;|^^. 

One of the particles in the perfect crystal, for example 
Qi, is set at the origin of the coordinate system. In 
general, p^'^\x, 0) depends on the choice of ai since each 
particle is not necessarily equivalent in the configuration 
of the target structure. In such a case, the average 
values of A over all possible choices of Oi should be 
used. However, this is not necessary for the target crystal 
structures used in this study. 

I have numerically determined A (a: 1 , 2 : 2 ) according to 
(9) using MG integration as implemented in the VEGAS 
algorithm in GNU Scientific Library. 
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FIG. 1. The density plot of A(xi,X2) for the square lattice. 
The location of the peak in A(xi,X 2 ) is first estimated on a 
coarse grid (Axi = 0.1) and is then estimated on a finer grid 
(Axi = 0.02) around the first estimated peak. 


TABLE I. The radial distances and coordination numbers for 
the first four nearest neighbors for perfect lattices. 


Square Honeycomb Kagome Triangular 


First 

1, 4 

1, 3 

1, 4 

1, 6 

Second 

\/2, 4 

\/3, 6 

Vs, 4 

Vs, 6 

Third 

2, 4 

2, 3 

2, 6 

2, 6 

Fourth 

^/5, 8 

\/7, 6 

V7, 8 

V7, 12 


The density plots of A(xi,X 2 ) for target lattices are 
shown in Figs. 1, 2, and 3. The peaks in A(xi,X 2 ) 
for square, honeycomb, and kagome lattices are located 
at (xi,X 2 ) = (2.00,1.04), (1.72,1.14), and (1.78,1.14), 
respectively. The associated LJG potentials are plotted 
in Fig. 4. All predicted LJG potentials have positive 
values at their first minimum; this reduces the number 
of the nearest neighbors and prevents the formation of 
an equilateral triangular lattice. The second minimum 
of the potential for the square lattice at a; « 2.19 
covers both the third and fourth nearest neighbors (see 
Table I for the neighbor distances and coordination 
numbers). Except around the second minimum, the 
potentials for the honeycomb and kagome lattices show 
similar behavior. The fact that these two lattices share 
the same second and third nearest neighbor distances 
likely contributes to the similarity in their potentials. It 
is reasonable that the second minimum of the honeycomb 
(kagome) lattice potential is located closer to the second 
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FIG. 2. The density plot of A(xi,X2) for the honeycomb 
lattice. The estimation grid is the same as in Fig. 1. 


(third) nearest neighbor distance, because the second 
minimum of the honeycomb (kagome) potential has a 
larger coordination number at the second (third) nearest 
neighbor distance than the third (second) one. 

The square lattice potential found in Ref. II has two 
minima at the first and fourth nearest neighbor distance, 
like the square lattice potential found here. But these 
potentials are not similar, because the former has a deep 
minimum at the first nearest neighbor distance and a 
very shallow minimum at the fourth nearest neighbor 
distance. The potential for honeycomb lattice here and 
that in Ref. 10 have some common features: the first 
minimum has positive value and the second minimum 
is at the second nearest neighbor distance. The potential 
for the kagome lattice here is completely different from 
the known potentials for the kagome lattice, which are 
purely repulsive. 

To examine the validity of these predicted potentials, 
I performed MG simulations for a constant number of 
particles N, volume (area) F, and temperature T. For 
each simulation, JV was greater than 500. The simulation 
box was a square of side L with periodic boundary 
conditions; L = ^/V was defined so that the density of 
the system N/V corresponds to that of the target perfect 
lattice. If N is not appropriate to fill the simulation 
box with the unit cells of the target lattice, the success 
rate of the crystallization into the target lattice will 
decrease. For square lattice potential, N was chosen to 
be a “magic number,” with which the unit cells fill the 
simulation box when one of the sides of the square unit 
cells and that of the simulation box are parallel. For the 
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FIG. 3. The density plot of A(xi,X 2 ) for the kagome lattice. 
The estimation grid is the same as in Fig. 1. 



FIG. 4. The predicted LJG potentials for square (dashed 
black curve), honeycomb (solid red curve), and kagome 
(dotted green curve) lattices. These are given by Eq. (10) 
with (xi,X 2 ) = (2.00,1.04), (1.72,1.14), and (1.78,1.14), 
respectively. 


honeycomb and kagome lattice potentials, TV was chosen 
so that the distortion of the lattice remains small when 
the target lattice fits in to the simulation box under the 
periodic boundary conditions. A random configuration 
was used as an initial configuration for each series of 
cooling process. Ten simulation runs were performed for 
square lattice potential, and all resulting solid phases are 
square lattice with some defects and no grain bound¬ 
aries. For the honeycomb and kagome lattice potentials, 



FIG. 5. Results of the 529-particle MC simulations for the 
LJG potential with (xi,X 2 ) = (2.00,1.04), annealed from 
ksT = 2.0 to ksT — 1.0 at N/V = 1 in the square simulation 
box of linear dimension L = 23.0. The particle pairs separated 
by a distance smaller than 1.2 are connected by a line segment 
to guide the eye. 

most simulation runs resulted in solid phases with grain 
boundaries and many defects; however, these solid phases 
exhibited local structures that resembled those of the 
target lattice. Several runs (six and four runs, out of 
thirty runs, for honeycomb and kagome lattice potential, 
respectively) resulted in crystalline structures without 
grain boundaries and only a few defects. The snapshots 
of these phases are shown in Figs. 5, 6, and 7, and the 
figures clearly show that the particles interacting with the 
LJG potentials predicted by the free-energy functional 
theory self-assemble into the target crystals. 

IV. DISCUSSION AND SUMMARY 

In this paper, I developed the free-energy functional 
method for the reconstruction of the interparticle inter¬ 
action potential from a given structure by combining 
the variational principle for the functional A[p^ "0] with 
Percus’s idea. In this free-energy functional method, the 
desired interaction potential is given as the function that 
maximizes "0]. This method was successfully 

applied to the square, honeycomb, and kagome lattices. 

The free-energy functional method introduced here 
requires single- and two-particle densities, p{x) and 
p^'^\x,x'), as input. To obtain the interaction potential 
corresponding to an artificially designed structure, e.g., 
the sets of atomic positions {fli}, p{x) and p^‘^\x,x') 
must be designed or approximated as was done in this 
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FIG. 6. Results of the 544-particle MC simulations for the 
LJG potential with (xi,X 2 ) = (1.72,1.14), annealed from 
ksT = 2.0 to ksT = 0.40 at N/V = 4/3\/^ in the square 
simulation box of linear dimension L — 26.58. 


FIG. 7. Results of the 504-particle MC simulations for the 
LJG potential with {xi,X 2 ) = (1.78,1.14), annealed from 
kBT = 2.0 to k^T = 0.4 at N/V — 3/2-v^ in the square 
simulation box of linear dimension L = 24.12. 


work. If, instead, the method is used to obtain a model 
potential that reproduces an experimentally observed 
structure, experimental data for p{x) and p^'^\x,x') are 
sufficient. 

The self-assembly of the target lattices in the MC 
simulations shows that the free-energy functional method 
works properly. However, this does not necessarily 
mean that the interparticle interactions obtained here 
are entirely optimized for the assembly of the target 
lattices. Indeed, if we perform several MC simulations 
using {x\,X 2 ) = (1.70,1.06), which is slightly different 
from the predicted parameters for the honeycomb lattice, 
the honeycomb lattice of comparable quality to the one 
shown in Fig. 6 is obtained more frequently than when 
the predicted parameter is used. This is likely due to 
the simplistic approximations employed here, such as the 
second order functional Taylor expansion in the activity 
for A in Eq. (9), the use of the Gaussian approximation 
for p{x) in Eq. (11), and the independent-fluctuation 
approximation for p^'^\x,x') in Eq. (12). More so¬ 
phisticated approximations will improve the resulting 
interparticle interaction potentials. The maximization of 
A in {xi,X 2 ) space may also give rise to the unsatisfactory 
quality of the resulting structures. The maximization on 
the full functional space of interparticle interactions is 
necessary for the best prediction within the framework 
of the present theory. 

The Taylor expansion of A not only affects the accu¬ 
racy of the prediction but also the range of applicability 
of this method. As discussed in Sec. Ill, the choice of the 
temperature and the trial potential is restricted to justify 


the truncation of the Taylor series. Therefore, the present 
form of this method does not provide a unified description 
of the optimal interaction potential for a given structure 
over a wide range (including T = 0) of temperature. 
Unfortunately, better approximations for Alp^ip] other 
than the Taylor series are not available yet. 

Even within the truncated Taylor series approxima¬ 
tion, the restriction on the temperature can be relaxed 
using non-negative trial potentials. Although this is an 
additional strong restriction on the interaction potential, 
we can expect that non-negative interaction potentials 
suffice for a wide variety of target structures at finite 
temperature, considering that such potentials produce 
various ground state structures. 

The interaction potential considered here is restricted 
to the isotropic one for simplicity. It was found to 
be sufficient for the self-assembly of the target crystals 
considered here, as expected from the past work for the 
ground state. The class of target crystals that cannot 
be assembled by isotropic particles is not yet clear,but 
the isotropy-assumption of the potential will break down 
and must be discarded if the target structure is strongly 
anisotropic. 

Improving the numerical efficiency is necessary for 
the future work, e.g., the use of the trial potential 
with many variational parameters; the inverse problem 
for the 3D structure, which requires time-consuming 
six-dimensional numerical integration in Eq. (9). The 
use of a terraced (discretized) interparticle potential is a 
possible candidate for the way to reduce the computation 
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time. The terraced potentials are used in molecular 
simiilations^®^^® and analytical calculations^®’^® as an 
efficient and satisfactorily accurate way to investigate the 
many body systems. 

The theory introduced here is also applicable to the 
inverse problem for liquids; in this case, the interparticle 
interaction is reconstructed from a given pair distribution 
function g{x,x') or a radial distribution function g{r). 
In fact, the method is more suitable for homogeneous 
simple liquids than for crystals, because for liquids, the 
approximations (11) and (12) for distribution functions 
are unnecessary. This is because p{x) = pis constant and 
p^‘^\x,x') is equal to p^g{x,x') in homogeneous simple 
liquids. While the inverse problem for liquids has been 
solved by the reverse MC method, our free-energy 
functional method serves as a new theoretical approach 
to tackle this problem. 
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